Propagation front of correlations in an interacting Bose gas 
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We analyze the quench dynamics ol a one-dimensional bosonic Mott insulator and focus on the 
time evolution of density correlations. For these we identify a pronounced propagation front, the 
velocity of which, once correctly extrapolated at large distances, can serve as a quantitative charac- 
teristic of the many-body Hamiltonian. In particular, the velocity allows distinguishing the weakly 
interacting regime, which is qualitatively well described by free bosons, from the strongly interacting 
one, in which pairs of distinct quasiparticles dominate the dynamics. In order to describe the latter 
case analytically, we introduce a general approximation to solve the Bose-Hubbard Hamiltonian 
based on the Jordan- Wigner fermionization of auxiliary particles. This approach can also be used 
to determine the ground-state properties. Complementary to the fermionization approach, we derive 
explicitly the time-dependent many-body state in the non-interacting limit and compare our results 
to numerical simulations in the whole range of interactions of the Bose-Hubbard model. 



I. INTRODUCTION 

In the past two decades, progresses in atomic physics, 
quantum optics and nanosciences have propelled quan- 
tum many-body theory to meet new challenges. It is 
indeed now possible to engineer systems that are con- 
crete realizations of some paradigmatic models, which 
were once introduced to grasp fundamental properties of 
more complex materials. New frontiers thus have to be 
explored, among which the dynamics of these isolated 
quantum models far from equilibrium is probably the 
least well understood and one of the most exciting. 

One of the fundamental questions that has to be ad- 
dressed is how correlations propagate in these systems. 
The Schrodinger equation allows in principle for corre- 
lations between distant points to build up in arbitrary 
short times [Tj. This contrasts with relativistic quan- 
tum field theories, where physical effects cannot propa- 
gate faster than the speed of light and causality relations 
between two points in space-time can only exist if one lies 
within the light-cone of the other. In a seminal work [2], 
Lieb and Robinson however showed that non-relativistic 
quantum many-body systems can still exhibit some sort 
of locality: In generic one-dimensional spin models with 
finite-range interactions, the propagation of correlations 
appears to be bounded by an effective light cone, out- 
side of which correlations are exponentially suppressed. 
Here the role of the speed of light is played by a velocity 
which is an intrinsic property of the many-body Hamil- 
tonian. The existence of so-called Lieb-Robinson bounds 
has many far-reaching implications. For example, they 
make it possible to simulate on classical computers the 
ground state properties as well as the dynamical evolu- 
tion of such quantum systems |3ji6j. They also provide 
a general link between the presence of a finite spectral 
gap and the existence of a finite correlation length in 
the ground state of certain lattice systems [7lI10|. How- 
ever, the extent to which Lieb-Robinson bounds can be 
generalized beyond spin systems remains an open ques- 



tion. Proofs or evidences for the existence of such bounds 
have indeed been reported in various systems, ranging 
from harmonic chains to the Bose-Hubbard model [lll - 
[T5] . But it is also possible to construct models in which 
the propagation velocity of correlations is explicitly un- 
bounded [IT] . 

Dynamical properties of correlations in a closed sys- 
tem can be probed by studying the time evolution fol- 
lowing a sudden change of parameter in the Hamilto- 
nian, a situation referred to as a quantum quench. The 
quench has a particular appeal in the context of ultra- 
cold gases in optical lattices as the relevant parameters 
in the Hamiltonian governing these systems can be easily 
varied in time |18j. Beside the existence of an effective 
light cone, it was discovered in recent theoretical studies 
that the time evolution of correlations in quenched sys- 
tems is characterized by a pronounced propagation front 
[n [la Hg Hg [ISHSS]. similarly to a Lieb-Robinson 
bound, the velocity at which this front propagates can in- 
volve a broad range of the spectrum of the Hamiltonian 
since the system is far from equilibrium. This makes 
the understanding of this feature particularly challeng- 
ing: Covariant low-energy effective theories would pro- 
vide a natural description [TTl [T^ EDH^ . but, due to 
the presence of high-energy excitations, realistic lattice 
models HH [El [El [ISl Hg [23H33] show a much richer be- 
haviour than their corresponding field theories. Gaining 
more insight into the non-equilibrium properties of quan- 
tum systems thus urges the development of new effective 
models. 

In a recent work |19| , the propagation of correlations in 
a quantum many-body system was studied both theoret- 
ically and experimentally in a one-dimensional bosonic 
gas in an optical lattice and a propagation front could 
be clearly identified. The observed behavior was inter- 
preted using an exactly solvable effective model derived 
from the Bose-Hubbard Hamiltonian and describing non- 
interacting fermionic quasiparticles. The key idea behind 
this model is to use a Jordan- Wigner transformation to 
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cure some of the problems inherent to the slave-boson 
methods proposed previously [34j [35] . In the present ar- 
ticle, we describe in more detail this new approach and 
use it to derive the ground state as well as the quench 
dynamics in the Mott-insulating phase. We show that 
its predictions are quantitatively correct in a regime of 
strong and intermediate interactions. Our model, being 
exactly solvable, allows us to explore the time evolution 
of the system at long times and we can show that the 
velocity of the propagation front exhibits a generic scal- 
ing behaviour. Using numerical simulations, we find that 
this behavior holds in all interaction regimes, down to 
the non-interacting limit of free bosons, where explicit 
solutions are available. The asymptotic value of the ve- 
locity of the propagation front, which strongly differs be- 
tween the strongly and the weakly interacting limits, can 
be used to characterize the crossover between these two 
regimes. 

This article is organized as follows: In Sec. |ll] we 
present the model that we will study, in Sec. |III| we carry 
out the fermionization procedure and derive general rela- 
tions that enable the calculation of equilibrium (Sec. IV) 



and non-equilibrium (Sec. |V]) properties. The velocity of 
the propagation front at weak and strong interactions is 



analyzed in Sec. |VIJ In Sec. |VII| we present our conclu- 
sions. 



II. ONE-DIMENSIONAL SYSTEM OF BOSONIC 
ATOMS IN AN OPTICAL LATTICE 

In this work we consider a one-dimensional system 
of bosonic atoms in an optical lattice. If the lattice is 
deep enough, this system can be described by the one- 
dimensional single-band Bose-Hubbard Hamiltonian: 



H 



^|-J(aJa,_ 



1 +h.c.) + ^(rij ~nf 



(1) 



where aj and aj represent the annihilation and creation 
operators of a bosonic atom at site j and rij 



counts the number of atoms at that site. We use a lat- 
tice constant aiat = 1 and the system is considered to 
be infinitely large and homogeneous. The kinetic part 
of the Hamiltonian is characterized by the hopping am- 
plitude J; The on-site interaction strength U is related 
to the s-wave scattering length. We work at fixed com- 
mensurate filling n, where the model exhibits a quantum 
phase transition between a superfluid phase at low in- 
teraction strengths U / J < {U / J)c and a Mott-insulating 
phase at large interaction strengths U/J > {U/J)c- At 
the specific filling n = 1, the critical value is given by 
(t//J)c 3.4 [36,1SZ1- The Bose-Hubbard model is non- 
intcgrable [351 [33] and exhibits complex many-body prop- 
erties; Especially, its non-equilibrium properties are far 
from being fully understood. 

In order to benchmark the analytical approaches, we 
will perform exact numerical simulations of model ([T]) 



by means of the density matrix renormalization group 
(DMRG) [iQl [41] , an algorithm based on matrix prod- 
uct states (MPS) |5]. While the DMRG algorithm 
gives highly accurate results for the ground state, time- 
evolution HMS] can only be calculated for relatively 
short periods of time. 



III. FERMIONIZATION APPROACH TO THE 
STUDY OF THE BOSE-HUBBARD MODEL 

In the following, we will describe in detail how the 
Bose-Hubbard model can be mapped onto an effective 
model of non-interacting auxiliary fermions. The proce- 
dure consists of four main steps: (i) The local Hilbert 
space is reduced to only three states and (ii) auxiliary 
bosonic operators are introduced that allow to switch be- 
tween these states (Sec. [HI A ); (iii) The auxiliary boson 
operators are fermionized by a Jordan- Wigner transfor- 
mation (Sec. IIIB); (iv) A constraint on the fermionic 



operators is relaxed so that the effective Hamiltonian be- 



comes quadratic and can be diagonalized (Sec. Ill C I 



A. Auxiliary boson representation 

In the Mott-insulating phase and away from the critical 
point, the local density fluctuations around the average 
filling n are limited. It is thus possible to truncate the 
local basis on a single site j to 3 states only: \n+m)j, with 
in = —1,0,1. Within this reduced basis, one can represent 
the bare atomic operators a^^^ in terms of constrained 

auxiliary boson operators 6^^]. with two flavors a = ±1 = 
±: 



(2) 



The '-I- '-bosons correspond to excess particles: b 



\n + l)^ , and ' — '-bosons to holes: &] 



\n)^ = \n-l)^. 



The local Fock state \n)j represents the vacuum state of 
the auxiliary particles bj,^\n)j = 0. Bosonic commuta- 
tion relations are obeyed: 



0, 



(3) 



which allow for the unphysical situation of single sites 
being occupied by more than one auxiliary boson. There- 
fore, the auxiliary operators have to fulfill the hardcore 
constraint 



(4) 



and double occupancies of different species need to be 
eliminated by imposing 



.6,,_=0. 



(5) 
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This representation in terms of doubly-flavored con- 
strained bosons is slightly different from the one used 
in slave-particle techniques [Ml [351 B5Hi5] , in which one 
introduces one auxiliary operator for each of the three 
local states \n + m)j and the number of auxiliary bosons 
per site is constrained to be exactly one. 



B. 



Fermionization 



The operator constraints Q and ^ are difficult to be 
ensured in general and one often resorts to mean field 
[3SH1HI and perturbative [Ml [351 HH] approximations. In 
the special one-dimensional case, it is however possible 
to use Jordan- Wigner transformations [49l [50] that allow 
on the one hand for the exact treatment of the hard-core 
constraint Q and on the other to suppress local pairing 
of auxiliary particles. 

Here, we follow the standard procedure of Jordan and 
Wigner [49] and introduce auxiliary fermion operators 
Cj^cr with number operators nj g- 
commutation relations 



Cj „Cj,cr and anti- 



Using non-local string operators, 



(6) 



Z. 



^3- — ^3- + ' 



(7) 



we relate the auxiliary bosonic operators to the fermionic 
ones: 



(8) 



The string operator Zj,„ counts the parity of the number 
of fermions accumulated over all sites j' < j (including 
the '-|- '-fermion on site if cr = — ) and obeys the relations 



Z 



(9) 



As a consequence, number operators within both 
fermionic and bosonic representations coincide: 

and the original atom number operator can be written as 

Hj — 71+ J- — j' + • (11) 

Due to the fermionic statistics, the hard-core condi- 
tions Q are satisfied automatically. The remaining con- 
straint ([5]) can be formally accounted for by the global 
projector 7-* — Y[j 'Pjj with Vj = (1 — 7ij_+7ij__) eliminat- 
ing states with both species on the same site. It is now 
possible to show that, within the truncated Hilbert space. 



the original Hamiltonian ([T]) can be exactly represented 
by the following fermionic model: 



j ^ 

- Jx/nin + 1) (4+4+1,- - Cj,+Cj+i,_) -^h.c. 

+ |(n,-++n,-_)|7'. (12) 

We note that the effective hopping amplitudes for the 
two different flavors differ by the bosonic enhancement 
factor of Eq. 



Exact diagonalization within the approximation 
of unconstrained fermions 



In practice, it is difficult to take care analytically of 
the projector "P. We will thus carry out the calculations 
in the approximation of unconstrained fermions (UF), 
"P — > 1, leading to a quadratic Hamiltonian that can be 
diagonalized exactly. We will see that the UF approxi- 
mation is justified because the main source of creation of 
double occupancies would be a local pairing mechanism, 
which, in the fermionic representation, is suppressed by 
the statistics of the auxiliary particles. 

The Hamiltonian ( 12 ) with 7-* = 1 can be rewritten in 
momentum space as 



-Ck,o 



J2 Mk)icl+cl,^-- c-k,-ck,+) , (13) 



with the bare dispersions 

E+{k) = -2J{n + I) cos(fc) + U/2 , 
(fc) = -2Jn cos(fc) + U/2 , 

and an antisymmetric pairing parameter 

A(fc) i2Jyyn{n + I) sin(fc) , 



(14a) 
(14b) 



(15) 



which obeys A(— fc) = — A(fc) — A*{k). In analogy to 
the Gutzwiller approximation [51H53j , the accuracy of the 
UF approximation can be estimated via the translation- 
invariant expectation value of the local projector 



P+- 



1 



{Pfit}) - (n,-+n,-_) 



(16) 



This quantity is a measure for the population of unphys- 
ical states and gives the order of magnitude of the er- 
ror in local observables (Due to translational invariance 
site indices of observables can be dropped). Addition- 
ally, we will study the quality of the relaxation of the 
constraint ^ by comparison to the numerically exact 
DMRG method. 
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FIG. 1. Quasiparticle dispersions (Eq. 21 and Eg. |22a[ ) at n = 1 for different interaction strengths (thin lines). Curvatures can 
deviate significantly from the cosine form of the strong coupling limit (thick lines). The width and the gap of the quasiparticle 
bands depend on the type of the quasiparticle. At U / J = 8, the gap of the '+'-particle closes, signaling the breakdown of the 
UF approximation. 



The quadratic Hamiltonian H^j-p can be diagonalized 
via a Bogolyubov transformation by introducing the 
quasiparticle operators 



7I - = - - v{k)ck.+ ■ 

The functions u{k) and v{k) fulfih the relations 

u{—k) = u{k) = u*{—k) , 
v{—k) = —v{k) = V* (k) , 

and are determined by the following expressions: 

-2iA(fc) 



u{k) 



v{k) 



cos at an 



1 + 



J' 



i sin ^atan 



E+{k) + E_{k) 



-2iA{k) 
E+{k)+E^{k) 



/2 



/2 



.2J^n(n+ 1) 



U 



sin(fc) + O 



(17a) 
(17b) 



(18a) 
(18b) 



(19a) 
(19b) 
(19c) 
(19d) 



We infer from the above equations that the '+ '-modes 
are excess particles each dressed with absent holes and 
the '—'-modes are holes dressed with absent excess par- 
ticles. This is particularly evident from the perturba- 
tive expressions (19b) and (19d|. We also note that the 



quasiparticle operators (17) can be interpreted as Dirac 
spinors [54] . 

Using the quasiparticle operators, the Hamiltonian can 
be written in the diagonal form 



(20) 



The dispersion relation of the individual quasiparticles is 
ea{k) ^ ~aJ coii{k) + hw{k) . (21) 



Here 2hLo{k) is the energy of a pair of two distinct types 
of quasiparticles with opposite momenta, which is given 
by 



2hw{k) = ^[E+{k)+E^{k)f +^/\{k)\ 



[/ - 2 J(2n + 1) cos(fc) + O 



(22a) 
(22b) 



The exact dispersion relations for different interaction 
strengths are displayed in Fig. [T] together with the first 
order expansion in J/[/ (strong coupling expansion) . One 
can see that the profiles rapidly differ from their limiting 
cosine shape as the interactions are lowered. Eqs. (14) 



and (22a I show that the width of the energy bands de- 



pends only on the hopping amplitude J and on the av- 
erage filling n (via the Bose enhancement factor), but 
does not depend on the interaction strength. At large in- 
teraction strengths, the energy gap is proportional to the 
interaction strength and the gap of the '-I- '-quasiparticles 
is strictly positive when the interaction is above a certain 
threshold: 



U/J>4{n + 1). 



(23) 



Below this threshold our UF approximation breaks down. 
For n = 1, the range of validity of our model is thus 
limited to U/J > 8, which is above the superfluid-to- 
Mott insulator transition {U/J)c ~ 3.4, but significantly 
lower than the mean-field transition {U/J)c ~ 12. A 
description of the phase transition might be achieved by 
introducing auxiliary operators on the basis of a coherent 
state representation (see e.g. but this goes beyond 

the scope of this work. 

The slope of the dispersion relations ^a{k) describes 
the group velocity of the quasiparticles. Of particular 
interest will be the relative velocity of pairs of quasipar- 
ticles of distinct types and opposite momenta: 



y{k) 



(24) 
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whose maximal value 



maxfe |v(fc)| 



(25) 



plays an important role in the characterization of the 
non-equilibrium properties. This maximal velocity corre- 
sponds to the point where the curvature of w(fc) changes 
sign. It is located at \k\ « 7r/2 at large interac- 
tion strengths and is shifted towards lower momenta at 
smaller interaction strengths, as can be seen in Fig[T] In 
the relevant interaction regime ( 23 1 , the maximum veloc- 
ity is well approximated by: 



^ 2J(2n+l) / _ 8n{n+l)J^ \ ( 



(26) 



In particular, one sees that Vmax is a decreasing function 
of U/J. 



For the strictly positive quasiparticle energies ( 23 1 , the 



ground state at a value of U and J is the quasiparticle 
vacuum 

\MU/J)) = n^'"'(fc)7fc,+7-fc.-l«) (27a) 

k 

= \{{u{k)+v{k)clJ_^^_)\n). (27b) 

k 

The ground state at infinitely strong interactions, i.e. the 
Fock state with n particles per site, represents the 
vacuum of the bare excess particles and holes. 

D. Local observables and correlation functions 

We summarize in this section some general properties 
of the correlation functions in the quasiparticle formalism 
that will be used later to derive ground state and non- 
equilibrium properties of the system. 

For the ground state (27), but also for the time- 



dependent wave functions ( 48 ) which will be introduced 



in section |Vj the only non- vanishing single-particle cor- 
relation functions are 



9d 



\'-j+d,a'^3,'^/ 



1 

2tt 



9d 



\C-j+d,CTCj,a) 
1 

2^ 



dfce *'"*(cfc,^c_fc,ff) , 



(28a) 



(28b) 



where a = —a and the thermodynamic limit has been 
taken. Possible time dependence (equal time) is implicit 
and expectation values are site-independent for the ho- 
mogeneous systems under consideration. We note that 
the correlations of the different types are equivalent: 



9d 



9d 



9d 



9d 



(29) 



Therefore, also the quasiparticle densities do not depend 
on the flavor and we can define a single density of exci- 
tations: 



= 2<?o' 



(30) 



Since the Hamiltonian is quadratic, correlations of the 
occupancies can be related to the single-particle correla- 
tions using Wick's theorem, which gives us 



f\ (7.(7 1 2 

'^^ \9d I 



{nj+d,a){nj^a') (31) 

(32) 



In the special case d = 0, the fermionic statistics, to- 
gether with the symmetry with respect to exchange of 
fermionic flavor, imply that the on-site correlator G^!°q = 

Iffd^ol^ vanishes and the local double occupancy factor- 
izes: 



(33) 



The density of excitations thus fully determines all local 
properties, including the atom number fluctuations: 



/ = ((%" - '^)^) = "-ex(l - f^ex/2) . 



(34) 



Atom correlations can be related to correlations of aux- 



iliary particles. Making use of Eq. (11), we can for ex- 



ample express atomic density correlations in the following 
way: 

Cd = {rijUj+d) - {nj){nj+d) (35) 
= E(Gr-Gr) (36) 



= -2(1.9^ 



(37) 



Rather than the density itself, ultracold atoms experi- 
ments with single-site resolved imaging |55l 156) can ac- 



cess the parity Sj 



The expression of parity 



correlations turns out to be similar to that of density 
correlations: 



Sd = (sjSj+d) - {sj){sj+d) 



(38) 
(39) 



Both density and parity correlations are particularly sim- 
ple to evaluate within the present approach. Correlations 
including the non-local string operator ([t]), such as the 
single-particle correlations (a]aj_|-d), can also be com- 
puted, but they require the evaluation of the Toeplitz 
determinant [57]. Interestingly, the fermionic string ([t]) 
is equivalent to the string operator recently measured by 
Endres et al. 1551. 



IV. EQUILIBRIUM PROPERTIES OF THE 
MOTT-INSULATING PHASE 

In this section we discuss the equilibrium properties 
of the Mott-insulating phase derived within the uncon- 
strained fermion approximation. 
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(a) atom mmibci' fluctuations (b) n.u. density coii'elations (c) density correlations 




U/J U/J 



FIG. 2. Ground state properties at n = 1. The numerical evaluation of the UF equations is compared with the strong couphng 
expansion, as well as with exact DMRG simulations of the Bose-Hubbard model with the local Hilbert space truncated at a 
maximum site occupancy of either A'^max = 2 or A^max = 6. (a) Atom number fluctuation / as a function of the final interaction 
strength U/J. (b) Nearest-neighbor density correlation Cd=i as a function of the final interaction strength U/J. (c) Density 
correlations as a function of the distance d. 



As argued in the preceding section, the observables are 
related to single-particle correlations ( 28 1 , which for the 



ground state (27) can be evaluated straightforwardly: 



number fluctuation / — nex(l ~ ?^ex/2) is evaluated nu- 
merically using ( 42 1 and compared to the strong coupling 
expansion (see also [55] ) 



{ckMC-k\s) = Sk,k'u{k)v{k) 



(40) 
(41) 



with the coefficients u{k) and v{k) given in Eq. ( 19 ). The 
local density of excitations pO| ) can thus be calculated 
from 



riex = / d/c v'^{k) . 

TT 



(42) 



In the case of strong interactions, one can also derive an 
explicit expression from the expansion (19b) and (19d) 
of the coefficients u{k) and v{k): 



2J2_ _ 
'^ex = T7^'^('^ + I) + O 



[/2 



{/4 



(43) 



Combining Eqs. (33) and (43) gives an estimate of the 



occupation of unphysical states (16 1, 



p+_-nL/4<(8(l + l/n))- 



(44) 



where the right-hand side stems from the evaluation of 
(43) at the lowest interaction considered (23). For n — 1, 



is less that 6% and we expect the error on local ex- 
pectation values to be of similar magnitude. With this at 
hand, we can now consider the behaviour of the density 
correlations in the ground state. In Fig. [2][a), the atom 



2P 



C/4 



(45) 



as well as to the results obtained from DMRG simulations 
with a truncation of the site occupancy to iVmax ~ 2 or 
-^max = 6 (system size is 256 sites, 400 DMRG-states 
are retained). The predictions of the UF approximation, 
both from the numerical integration of (42) and from 



the strong coupling expansion, are in excellent agreement 
with the DMRG simulations for all interaction strengths 
satisfying ( 23 1 . The accuracy of the truncation of the 
local Hilbert space to 3 states only is also confirmed by 
the DMRG simulations. Higher occupancies start to be 
important only for interaction strengths below the point 
where the UF approximation breaks down. 

We further compare our results with the ones derived 
within a Holstein-Primakov approximation of the slave- 
boson representation used, e.g., by Huber et al. [5S]. This 
approach is equivalent to the auxiliary boson representa- 
tion ([2| when fully relaxing the constraints (4[5). We find 
that the local observables cannot be well described at in- 
termediate interaction strengths, even though the density 
of excitations is small (Fig. [2]) . A similar instability has 
been observed with slave bosons in Ref. [47^. We will 
analyse the slave-boson approach in more detail later in 
the context of the non-equilibrium dynamics (Sec. VD|. 

We can also evaluate analytically non-local density cor- 
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relations to second order in J/U: 

J2 



C/2 



J4 

C/4 



(46) 



As shown in Fig.|2jb), the above expression only slightly 
overestimates the amplitude of the correlations compared 
to the full numerical evaluation of the integrals ( 28 ) and 
(37) with (41 1. We therefore conclude that local ob- 



servables and nearest-neighbor correlations in the Mott- 
insulator regime are well described by a perturbation ex- 
pansion to order ,P /U^ . This is no longer the case for 
longer-range correlations, which are simply vanishing ac- 
cording to the expansion to second order, whereas the 
exact DMRG predicts that they should be finite and ex- 
ponentially decaying. As can be seen in Fig. [2f^c), the nu- 
merical evaluation of the UF equations provides a much 
better agreement with the DMRG results. The correla- 
tions at d = 2 can be almost perfectly reproduced and 
a similar decay length is found. The amplitude of the 
correlations for d > 2 is overestimated, however, and the 
discrepancy becomes worse as d increases. 



For the wave function ( 48 1 , the non- vanishing equal-time 



single-particle correlations evaluate to 

+Sk,k'v{k)v{k) \e^"^^'''>'v{k)u{k) - v{k)u{k) 



(50a) 



and 



{clA*)ck'.a{t)) = Sk,k' [2cos{2uj{k)t)u{k)u{k)v{k)v{k) 

u^{k)v^{k)-v'^{k)u^{k)\ . 

(50b) 

Based on these expressions, the expectation values of any 
observables can either be calculated analytically, when 
the strong coupling expansion holds, or computed nu- 
merically for lower interactions (see III D ) . 



B. Strong coupling expansion 



V. QUENCH DYNAMICS - GENERAL 
DESCRIPTION 

We analyse the quench dynamics of a system prepared 
initially in a deep Mott-insulating state. We first derive 
the general results for the time evolution of the wave 
function and the correlation functions and then give ex- 
plicit expressions for the case where the initial state is a 
Fock state (infinite interactions). These results form the 
basis for the detailed discussion of the physical properties 
of the quench dynamics in the subsequent section [yi] 



A. Time-dependent wave function and correlations 

The initial state considered is the ground state at some 
values of J and U satisfying the condition ( 23 1 and takes 
the form 

iV-in.) = n ("o(^) + ^o(fc)cl.^+cl, _) \n) . (47) 

k 

The time-evolution of this state under the Hamiltonian 
i?UF reads 



-iHupt/h 



with 



n {m - me~'"'^'H,+i-k,-) \Mu/j)) , 

k 

(48) 



u{k) = u{k)uo{k) - v{k)vQ{k) , (49a) 
v{k) = vik)uoik) - uik)vo{k) . (49b) 



For concreteness, we focus now on a quantum quench 
starting from the Fock state with filling n by setting 
wo(fc) — 1, vo{k) = and thus u{k) = u{k), v(k) = v{k). 
In this case, the expansion in J/U leads to 



/ / \ I \\ U Jn{n -f 1) . 
{Ck,a[t)C-k,s{t)) = I sm(fc) 



-2iu}{k)t 



1 



(51a) 



{cljt)ckAt)) = 



8J^n(n + l) 
C72 



sin^(A;) [cos(2a;(fc)i) - 1] 



(51b) 



where aj(fc) stands for the dispersion in the strong cou- 
pling expansion ( 22b ) . 

Within this expansion, the dynamics of the lo- 
cal density of excitations can be expressed in terms 
of the Bessel functions of the first kind, Jn{z) ~ 
r dfce-*^'=°''('=)+"'=. One gets 



"-ex(i) 



8fi(n-h 1)J2 



t/2 



1 — cos 



{ut/h) (j2{Jt) + MJt) 

(52) 

with J — 2J(2n -|- l)/h. In the relevant interaction 

regime, the population of unphysical states p_| (t) thus 

remains as small as in the ground state ( 44 1 and we ex- 



pect the UF approximation to be well behaved in general. 
It is, however, important to note, that the expansion is 
not rigorous since the approximative dispersion ( 22b ) is 
multiplied by time, which is unbounded. 
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The single-particle correlators required to derive non- 
local density correlations (37) read 



5r(t)-^'^'^f|^/rffce'-sin(fc) 



'2iLj(k)t 



1 



+Jd-i{Jt)) + Sd 



Making use of the identity 



2d 



Jd+iiz) + Jd-iiz) = —Jdiz) , 



1 • 

(53) 



(54) 



we obtain the following expressions for the non-local den- 
sity correlations: 



(55a) 
(55b) 



Cd>i{t) 



We note that the interaction strength U is involved 
only in the magnitude of the correlations for d > 1, via 
the dimensionless parameter J/U. In the case of nearest- 
neighbor correlations, we find an additional oscillation of 
the amplitude with frequency U/h. 



C. Accuracy of the UF approximation 

In this section, we analyze the accuracy of the succes- 
sive approximations that lead from the Bose-Hubbard 
model to the UF approximation and its strong coupling 
expansion. For this purpose, we introduce the root-mean- 
square differences 



Xd = 



(t) " 



(56) 



where C^^'' and C^^-* are the density correlations pre- 
dicted using two different level of approximations. By 
observing the dependency of Xd with the distance c?, we 
can verify, in particular, whether a given approxima- 
tion breaks down at large times. We use imax = 3/i/J, 
the maximal time accessible by our DMRG simulations. 
We use a DMRG algorithm in the thermodynamic limit 
[BOt l6T] . with a second-order Suzuki-Trotter decompo- 
sition of time step At — Q.02h/U, and we retain 2400 
states. The numerical error is always smaller than sym- 
bol size and line width. 

In Fig . [3[ a), we first compare the strong coupling ex- 
pansion ( |55p and the numerical evaluation of the UF ap- 
proximation. We find that the expansion is relatively 



accurate {Xd < 10""^) down to interactions U/J ~ 10, 
except for the d = 2 correlation, which only slowly con- 
verges to the exact results when increasing U/J. We 
remind here that a similar accuracy is reached for the 
ground state correlations (Fig.[2|. 

In Fig.|3][b), we then compare the numerical evaluation 
of the UF approximation with the exact DMRG simu- 
lation of the Hamiltonian (12), i.e. the Bose-Hubbard 



model when the site occupancy is truncated to iV,nax — 2 
The UF approximation appears to be accurate within 
xl < 10"^ for U/J > 12. As we will show in Sec. |VIC 
the UF approximation still qualitatively describes the dy- 
namics between U/J~ 12 and U/J = 8, which marks 
the break down of the quasiparticle picture. 

Finally, we compare in Fig. [sj^c) the predictions of the 
DMRG simulation when the local Hilbert space is trun- 
cated to a maximum site occupancy A'max = 2, corre- 
sponding to the model (fT2|), or A'^,„ax = 6. We observe 



that the error due to the truncation starts to be sig- 
nificant {Xd > 10~^) only for U/J < 6, i.e. when the 
interaction energy becomes larger than the width of the 
quasiparticle band. 



D. Comparison with the Holstein— Primakov 
approximation for auxiliary bosons 

In order to generalize the description to higher dimen- 
sions, it may appear tempting to fully relax the hardcore 
constraint ^ and work with bosons instead of fermions. 
The resulting Hamiltonian is then equivalent to the one 
derived by Huber et at [35] using Holstein-Primakov 
bosons. One then obtains equations very similar to those 
derived in the fermionic model, except that the coefficient 
v{k) becomes symmetric instead of antisymmetric. As a 
consequence, local pairing of different species is no longer 
suppressed and the occupation of the unphysical states 
become much larger that in the fermionic approach. In 
order to quantify this effect, one first has to calculate the 
single-particle correlations. To lowest order in J/U, we 
find 



j cos{k)t 



(57) 

which can be compared to the fermionic version ( 53 1 . It 



follows from this equation that the overcompleteness 



2n(n-|- 1)J^ 
IP 



Ji{Jt) 



(58) 



is of the same order as the density of excitations and the 
density correlations. The Holstein-Primakov approxima- 
tion therefore fails to describe the quench dynamics, at 
a qualitative level even in the limit J7 3> J. This can 
be observed, for example, in the density correlation func- 
tions, that behave in a radically different way than their 
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FIG. 3. Root-mean-square differences of the density correlations (56 1 obtained from different degrees of approximation, 
(a) Tlie strong coupling expansion is compared to the numerical evaluation of the UF approximation, (b) The numerical 
evaluation of the UF approximation is compared to the exact DMRG simulation of the Bose-Ifubbard model with a local site 
occupancy truncated at A'^max = 2. (c) The DMRG simulation with A'^max = 2 is compared to the DMRG simulation truncated 

at A^max = 6. 



fermionic counterparts: 



we can derive explicit equations for the density correla- 
tions: 



Cd>i{t) 



HP f n{n + 1) J 



(jd~i{Jt) - Jd+i{Jt)y 



(59) 



E. Limit of non-interacting bosons 



In this section, we complement the preceding analy- 
sis of quenches on the strongly interacting side by the 
extreme situation of a quench from infinite to zero in- 
teractions [33l [62] ■ At U/J = 0, the time evolution is 
readily described in the Heisenberg picture, 



in which individual bosons propagate with free disper- 
sion. In the thermodynamic limit, this yields the propa- 
gator 



2tt 



dkexp [—i {2J cos{k)t/h — kd)] 
2Jt 



(60) 



Using the following relation for the initial Fock state: 



{n\alaqalas\n) = n^Sp^qSr,s + n(n -I- 1)(1 - Sp^q)6p^sSq,r ■ 



Cd{t)= n2(^J^2(2ji/;,)j 

J 

+ fi{n + l)(^J2j,+di2Jt/h)Jji2Jt/h)] 

j 

- n{n +l)Y^ jf+^{2Jt/h)jf{2Jt/h) - 

3 

= -n(fi + l)Y,JfW^Jtlh)jf{2Jt/h) . (61) 



Here we have used the properties Jn{u ± w) = 
E™=-oo Jn^MJnM and JM = for d ^ 0. 



VI. HOW QUASIPARTICLE PAIRS CARRY 
DENSITY CORRELATIONS ACROSS THE 
SYSTEM 

We now analyse in detail how correlations spread in 
the quench dynamics starting from the Fock state with 
n atoms per site within the Bose-Hubbard Hamiltonian 
([T]). The description in terms of fermionic quasiparticles 
for intermediate and strong interactions ( 48 1 , as well as 



the non- interacting solution (61 ), provide a firm basis for 



the interpretation of the outcome of the DMRG simula- 
tions and of recent experimental results pjj and allow 
for their extrapolation at long times, where no analytical 
solution is available so far. 



A. Quasiparticle pairs 

For concreteness, we restrict our discussion in the fol- 
lowing to the filling n = 1, where the '-t- '-quasiparticles 
(17a I correspond to doublons and ' — '-quasiparticles 
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(17b) to holons. The relevant processes involved in the 



quench dynamics can be best understood in the expan- 
sion of the wave function ( 48 ) to lowest order in the aux- 
iliary fermion operators: 

IV' W) = In) +^^5]sin(fc)4,^^,^Jf^) (62a) 

k 

z^5:sin(fc)e^«'^™^('=)*/'-'ct_^cL,,_|n). 



U 



(62b) 



In this representation, the state decomposes into two 
parts: a time-independent part (62a) consisting of 
the Fock state and the symmetric superposition of 
bound nearest-neighbor doublon-holon pairs, and a time- 
dependent part (62b) describing the superposition of 
propagating doublon-holon pairs. The dynamics fol- 
lowing the quench is driven by the propagating pairs, 
whereas the steady state is solely determined by the 
bound pairs, as the contribution of the propagating pairs 
phases out at long times. At lowest order in J/J7, the 
steady state simply corresponds to the ground state at 
the final interaction strength. Higher order terms in the 
strong coupling expansion would describe the population 
in the excited states. At t = 0, the bound and propa- 
gating pairs interfere destructively and one recovers the 
initial Fock state \n). Finally, we note that the wave 



function ( 62 ) is equivalent to the one obtained within a 
time-dependent perturbation theory in Appendix [A] and 
can be used to derive the perturbative results presented 
in Sec. IVBl 

The doublon and the holon forming a propagating pair 
are produced initially on neighboring sites by a single 
hopping event and then move in opposite directions. The 
two quasiparticles are entangled, since the pair is de- 
scribed by a superposition state: 



This ensures a constant atomic density and leads to 
strong bipartite entanglement as the pairs are stretched 
across the system |TT]. The momentum distribution of 
the quasiparticle pairs is sine shaped, from which follows 
that the quasiparticles propagate as a wave packet. The 
maximal weight of the momentum distribution is located 
at the wave vector |fc| = 7r/2, where the dispersion rela- 
tion (22b) is close to being linear, and is characterized 
by the maximal group velocity Vmax = 6J/h. The wave- 
packet structure of the propagation can also be made 
explicit by turning the sum over the momenta in ( 62b ) 
into a sum over the lattice sites: 



i6J cos(k)t/h J t 



\n) 



In the above expression, one immediately recognizes the 
propagation velocity 6 J/ h in the argument of the Bessel 
functions. However, we expect a large dispersion of the 
wave packet due to the width of the momentum distri- 
bution. A detailed description of the propagation of the 
quasiparticle pairs is left for Sec. |VI C| 

The situation is somewhat different for weakly inter- 
acting bosons. In the non-interacting solution ( |61[ ), the 
correlation functions result from the interference between 
free bosons propagating with a relative velocity AJ/h. 
Unlike the auxiliary particles in the strongly interacting 
case, the number of free bosons per site is not limited, 
which leads to some qualitative differences that we will 
discuss later. 



B. Correlation signal in the density correlations 

The equal-time density correlations Cd(t) in the 
strongly interacting limit exhibit a very peculiar feature 
for d > 2, namely the presence of negative signal, a dip, 
propagating to larger distances at longer times. This can 
be seen for example in Fig. [4j where we display the time 
evolution of these correlations for U/J~ 18, as predicted 
by the UF approximation and the DMRG simulation 
(which are in remarkable agreement). This character- 
istic signal is already present in the perturbative result 
and can be attributed to the propagating quasiparticle 
pairs, illustrating the interest of this picture. 

The structure of the nearest-neighbor correlation is 
more complicated. In the long-time limit and within the 
strong coupling expansion, the nearest-neighbor correla- 
tion reaches the value corresponding to the ground state 
at the final interaction strength. At short time, it ex- 
hibits oscillations driven by the interaction strength U 
and corresponding to the interaction of a holon (dou- 
blon) of the bound pair (62a) with a doublon (holon) of 
a propagating pair ( 62b ) . In the first order of the strong 
coupling expansion (62), the bound pairs extend only 



E 



M6Jt/h)cl^c]^, _\n) . (63) 



3Jt 



over a distance d = 1 but in the full numerical integra- 
tion they can spread over larger distances [cf. Fig. [2j(c)], 
leading to additional oscillations for d = 2. These oscilla- 
tions are clearly visible when evaluating numerically the 
correlations within the UF approximation, as well as in 
the DMRG simulations (Fig.|4|. 

While the UF approximation is almost exact at short 
distances, it overestimates the weak oscillations with pe- 
riod h/U at larger distances. We found that these stem 
from terms of order J'^/U^ dominating the doublon- 
doublon and holon-holon correlations. These oscillations 
are also present in the DMRG simulations, but with a 
much lower amplitude. Interestingly, these terms cancel 
in the parity correlations studied in 19J, where the UF 
approximation is in even better agreement with the exact 
simulations. 

For quenches to intermediate values of the interaction 
strength, the dynamics of density correlations exhibits es- 
sentially the same features as described above. This can 
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be seen in Fig. [5j where the dynamics fohowing a quench 
to C//J = 9 is depicted. In particular, the characteris- 
tic dip corresponding to the propagating quasiparticles is 
still present. We note that the propagation of this corre- 
lation signal is still remarkably well described by the UF 
approximation, even though this model is close to break- 
ing down at this interaction strength. However, one sees 
that the strong coupling expansion significantly overes- 
timates the amplitude of the correlations at d = 1 and 
that the amplitude of the unphysical oscillations in the 
numerical evaluation of the UF approximation increases. 

In the weakly interacting regime, one could expect a 
different behavior since the relevant quasiparticles are of 
different nature. However, the main features that char- 
acterize the dynamics of density correlations at strong 
and intermediate interactions are remarkably preserved, 
as can be seen in Fig. |6] In particular, a propagating 
dip can be identified in all cases. The main difference 



between the non-interacting (61) and the strongly in- 



teracting ( 55 ) cases is the lower velocity and very slow 
decay of correlations at long times when U = 0. At 
U/J — 2, this long tail is already strongly suppressed 
at short distances, but is still visible at longer distances. 
For U/J = 4 (Fig. [6]) , one sees that the overall profile of 
the propagating correlation signal is already very similar 
to the more strongly interacting case (Figs. [4]and[5|. 

At low interaction strengths, one has to be careful 
when using the DMRG simulations since the truncation 
of the local Hilbert space to a finite number of states can 
introduce significant errors. By comparing DMRG simu- 
lations to the exact formula (61 1 obtained in the "worst" 
case U/J = 0, we found that a maximal site occupancy 
A'max — 6 represents a fairly safe approximation, whose 
accuracy improves with the distance for the times con- 
sidered (see Fig. [6]). 

To summarize the analysis conducted in this section, 
we observe that the dynamics of the density correlations 
is dominated by the propagation of a negative signal 
(dip). For strong interactions, this dip results from the 
propagation of the quasiparticle pairs described in the 
section [yl} In the next section we will investigate in more 
detail the shape of that signal and characterize its prop- 
agation velocity quantitatively. 



C. Analysis of the signal propagation 




FIG. 4. Dynamics of density correlations at distance d > 1 
after a quench from the Fock state |n) at infinite interactions 
to a final interaction U/J = 18. The results for difi'erent dis- 
tances d are shifted for clarity by 0.005(d — 1). We display 
the results obtained from the numerical evaluation of the UF 
equations, from their strong coupling expansion and from ex- 
act DMRG simulations. The shaded blue profiles figure a 
Gaussian fit of the correlation signal from the DMRG simu- 
lation, after the high-frequency oscillations have been filtered 
out. The filled blue circles mark the center of the fitted pro- 
file, i.e. the position of the signal. The filled red circles mark 
the position of the correlation signal obtained in the same way 
from the numerical evaluation of the UF approximation (fit 
not shown). The Airy function that appears in the analyti- 
cal formulas derived from the UF approximation is plotted as 
inset. 



We can get a lot of insight into the propagation of the 
signal from the following approximation of the density 
correlations at large distances (we recall that the lattice 
constant aiat is set to one): 



Ai 



i^ (^-{2/dy/^6Jt/h~d)'^ 



(64) 

In the above expression, derived from ( 55b I , we made use 



of the relation existing between the Airy function Ai(— z) 



and the high-order Bessel functions [55] : 

Jd{d + zd^/^) - 2i/3d-i/3Ai( - + O(d-i) . (65) 

The Airy function is plotted in the inset of Fig. |4j It 
exhibits a peak located at zq ~ 1.02 and surrounded 
by an exponential tail on the side z < zq and by an 
algebraically-decaying oscillating tail on the side z> z^. 

Disregarding the monotonously and slowly varying 
prefactor in Eq. (64), the profile of the Airy function 
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0.30 




FIG. 5. Dynamics of density correlations after a quench from 
the Fock state |n) at infinite interactions to a final interaction 
U/J = 9. See Fig. [i] for more information. 



alone allows us to understand several features of the prop- 
agation of the correlation signal. For example, it reveals 
the existence of a well defined propagation front, since 
the correlations are exponentially suppressed for times 
t < ipcak, with 



J ^pcak 



1/3 



(66) 



The signal in the density correlations corresponds to the 
peak of the Airy function. Once this peak has passed, 
that is for t > ipcakj the correlations show an algebraic 
decay with oscillations. From the definition of ipoak, 
one sees that two terms contribute to the propagation 
of the correlation signal: the first is simply proportional 
to the distance, corresponding to a well defined velocity, 
whereas the second is proportional to d^^^. The linear 
contribution dominates at large distances, leading to a 
light-cone-like spreading of the correlations. At small 
distances, however, the dynamics deviates significantly 
from the asymptotic light cone. This behavior can be 



FIG. 6. Dynamics of density correlations after a quench 
from a Fock state |n) at infinite interactions to weak final 
interactions. The results for different distances d are shifted 
for clarity by 0.25(d— 1). Unlike in Figs. |4] and [5] the position 
of the correlation signal of the DMRG results is identified 
with the position of the absolute minimum and it is denoted 
by the circles in the corresponding colors. 



accounted for by defining an 'instantaneous' propagation 
velocity: 



Vd = [tpcak{d + 1) - Vak(rf)] ^ 

= v^(l-^d-^/3)+0(d-^/3). (67) 



One sees immediately in the above equation that the 
asymptotic light cone is characterized by the velocity 
Voo = QJ/h and is reached algebraically at large dis- 
tances, whereas the propagation velocity can go down 
to approximately 4:J/h at short distances. 

A similar analysis can be carried out for the case 
U/J = 0. It turns out that the correlation dip is almost 
completely described by a single term in the infinite sum 
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FIG. 7. (a) Instantaneous propagation velocity in the UF approximation or in the non-interacting case as a function of the 
distance d (filled symbols). Light symbols represent DMRG data. Lines show the fits jv^ — Voo| oc The data have been 

shifted vertically for a better visibility, (b) Instantaneous propagation velocity obtained by DMRG simulation as a function of 
the distance d (filled symbols). Lines show the fits |vd — Vooj oc d~" with fixed exponent a — 0.65. (c) Asymptotic velocities 
Voo extracted from the finite distance data versus interaction strength using \vd — Voo| oc d""'^^. Error bars denote the 2-sigma 
uncertainty of the fit that yields the asymptotic velocity. 



(68) 



(61 1. For even d, for example, we obtain: 

= -2d-^/^M'^ {Ut/h - d) /2) 



The same expression for the instantaneous velocity ( |67[ ) 
therefore holds in the non-interacting case as well, but 
with Voo = 4 J/fi,, which is the velocity of freely propa- 
gating bosons. The behaviour at U/J = mostly differs 
from the strongly interacting case once the correlation 
dip has passed {t > tpcak): further terms (61) beyond 
Eq. (68) then become important, which causes correla- 
tions to decay very slowly (see Fig. |6]) . 



The width and the height of the correlation dip can also 
be derived from the expressions ( |64|68[ ). For both the 
interacting and the non-interacting case, the width in- 
creases proportionally to d^^^ while the height decreases 
with d~^/^ /U"^ in the strongly interacting case and 
with d~^^^ for U/J = 0. We note that similar power 
laws have been found for the quantum Ising model |23) . 

In the following, we show that the approximate scaling 
of the velocity found in the strongly and non-interacting 
limits holds for any interaction strength. We first con- 
centrate on large interaction strengths. Within the UF 
approximation, we can evaluate the correlations up to ar- 
bitrarily long times and make a rigorous scaling analysis 
of the instantaneous propagation velocity. We determine 
the position of the dip by means of a Gaussian fit af- 
ter having filtered out oscillations with a period shorter 
than h/U using a low-pass filter. Fig. [7][a) illustrates 
for a few interaction strengths that the analytical scaling 



behavior \vd — ^oo\ oc d~°' is accurately reproduced at 
sufficiently large distances c? > 5. Extracting the asymp- 
totic velocities Vqo and the exponents a with a fit over 
distances 6 < d < 400, we obtain values in very good 
agreement with the approximated analytical predictions. 
For example, the exponent is found to be the same for 
all interactions: a = 0.650 ± 0.002. The small difference 
to the value a = 2/3 expected from the Airy functions 



(67) is most probably due to the prefactor in (64), which 



we neglected when deriving (67). The asymptotic veloci- 



ties match the ones that we expect from the quasiparticle 
dispersion relation Eq. ( 25 ) , as shown in Fig.jTjc). Close 
to the breakdown of the UF approximation, the oscil- 
lation frequencies due to the interaction and the finite 
bandwidth become similar and one cannot easily filter 
out the first one anymore. The instantaneous velocity 
Yd therefore shows an oscillatory behaviour even at very 
large distances d < 50. Nevertheless, the scaling be- 
haviour remains perfectly obeyed on average and in the 
long-distance limit. In the non-interacting case, shown 
in Fig. [Tjja), we extract accurately the position of the 
correlation signal by simply locating the first minimum. 
We again find the scaling exponent a — 0.650±0.002 and 
the extracted asymptotic velocity is close to the expected 
value Voo = ^J/h [cf. Fig.[7];c)]. 

Using DMRG simulations, we can calculate the dy- 
namics exactly for all interaction strengths, but we are 
restricted to short time and length scales. We therefore 
fix the scaling exponent to a = 0.650 in order to ex- 
tract the asymptotic velocities. In Fig.[7][b) we show that 
the scaling |vd — Voo| oc d~^-^^ becomes accurate as the 
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distance increases for both strong {U / J > 8, extracted 
with low-pass fiUer and Gaussian fit) and weak interac- 
tions {U/J < 4 extracted directly from the peak with- 
out low-pass filter). Despite the limited number of data 
points available in the scaling region, we can determine 
the asymptotic velocities with a reasonably small uncer- 
tainty. The values that we obtain, gathered in Fig. T^c), 
are in good agreement with those predicted by the UF ap- 
proximation. The lack of data in the range 4 < [// J < 8 
results from the mixing of the time scales related to ki- 
netic and interaction processes and which prevents us 
from locating accurately the position of the correlation 
signal. The asymptotic velocities in Fig.jTf^c) can be seen 
as a characterization of a crossover between a regime of 
quasi-free bosons {U/J < 4), with a renormalized veloc- 
ity, and the strongly interacting regime described by two 
flavors of fermions. This crossover is not directly related 
to the ground-state phase diagram of the Bose-Hubbard 
model since the propagation velocity reflects the disper- 
sion in the center of the Brillouin zone (at wave vectors 
k « ±f ), rather than low- wavelength modes. As a conse- 
quence, Voo is considerable higher than the sound velocity 
in the superfluid regime |641 165| and a linear propagation 
with Voo ^ 6J/h is found at strong interactions, even 
though at equilibrium the system would be in the Mott- 
insulating phase. 

As a final remark, we note that the dependency of the 
propagation velocity on U/J in that system has been 
studied before by Lauchli and KoUath [15^, who con- 
sidered the case of a quench from a small interaction 
strength to a larger one. Surprisingly, the instantaneous 
spreading velocity has been found to exhibit a maximum 
at intermediate interaction strength. A possible expla- 
nation for this effect could be that bosonic atom number 
fluctuations present in the initial superfluid state may 
lead to enhanced velocites as compared to the quench 
from the Fock state. A quantitative comparison between 
our predictions and this previous work would require an 
extrapolation of the velocity to large distances which is 
difficult in the absence of an analytical model. 



VII. CONCLUSIONS 

In order to describe accurately the quench dynamics of 
the one-dimensional Bose-Hubbard model in the Mott- 
insulating regime, we have developed a new analytical ap- 
proach relying on the fermionization of auxiliary bosons. 
Its predictions regarding both the ground state and the 
dynamical properties are found in quantitative agreement 
with exact numerical simulations for large and interme- 
diate interaction strengths U/J>8. This constitutes a 
great improvement with respect to the analytical models 
introduced previously. 

Using this model, we are able to investigate the time 
evolution of density correlations in the quenched system 
over exceedingly long times. We observe a characteristic 
light cone dynamics, meaning that there exists a distance. 



linearly growing in time, beyond which correlations be- 
tween distant sites are exponentially suppressed. More 
precisely, correlations spread as a wave packet along this 
light cone, forming a propagation front whose position 
can be unambiguously identified. A careful analysis of 
the velocity with which this front propagates reveals a 
generic scaling behavior characterized by a universal ex- 
ponent and an asymptotic velocity dependent on the in- 
teraction strength. The same behavior is found in the 
non-interacting limit of freely propagating bosons, where 
an explicit solution is available, as well as in the inter- 
mediate regime < U/J < 8, where we rely on exact 
numerical simulations. The asymptotic velocity, which 
varies significantly between the weakly and the strongly 
interacting regime, is a useful quantity to characterize a 
broad spectral range of the Hamiltonian as it does not 
depend only on its low-lying modes. 

Building upon this first success, we envisage that the 
representation of the Bose-Hubbard model in terms of 
fermionic quasiparticles could shed new light into the 
mechanism for thermalization or serve as a tool to in- 
terpret the outcome of spectroscopic measurement on 
laboratory systems, such as modulation spectroscopy or 
Bragg spectroscopy for ultracold gases in optical lattices. 
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Appendix A: Perturbation theory 

In this section we develop a complementary pertur- 
bative approach to recover the behaviour in the strong 
coupling limit to first non- vanishing order in J/U . The 
situation we consider is the quench from the Fock state 
\ip{0)) = \n) at filling n to a large final interaction 
strength U/J. 

In all generality, the wave function after a sudden 
change of parameters can be written in the eigenbasis 
\(j)n) (with corresponding eigenenergies En) of the final 
Hamiltonian 



-(0„|i^(O))|0„ 



(Al) 



Usually the difficulty lies in determining the eigenstates 
\<pn) and their corresponding energies i?„ in a many-body 
problem. In this appendix we determine \(pn) and En by 
perturbation theory in J/U in a system of length L with 
periodic boundary conditions. Note, that this is not a 
full perturbative expansion, since we will not expand the 
exponential in the corresponding power. 
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We consider the interaction term of the Bose-Hubbard 
Hamiltonian as the unperturbed Hamihonian and the 
kinetic part as the perturbation. The eigenenergies of 
the unperturbed Hamihonian are muhiples of the in- 
teraction strength U and the corresponding states are 
Fock states. More precisely, the ground state is the 
Fock state \n) with vanishing energy. The lowest ex- 
cited states are the states with a single particle-hole ex- 
citation with energy U. These we denote by \(j){m,d)) 
with an occupation n for all the sites except for site m 
with n + 1 atoms and the site m + d with n — 1 atoms, 



-ibl 



Using degenerate 



V "("+!) 

perturbation theory (restricted to the same symmetry 
sector as the initial state) at first order in J/U, the 
ground state energy remains zero and the ground state 
of the final Hamiltonian is given by 

|0o) « \n) + + J2 i\Hm,l)) + |0(m, - 1))) . 



U 



m=0 



(A2) 



The lowest band of excited states resulting from the sin- 
gle particle-hole excitations is formed by 



^2n(n+ 1)J 



Here 10°) (p = 0, . . . ,L — 1) are the symmetric states that 
diagonalize the kinetic part of the Hamiltonian given by 



V2 
L 



L-l L-1 



sm{T:pd/L)\<j){m,d)) 



d=l m=0 



Note that the index d only starts at 1 to avoid the dou- 
ble counting of the Fock state. We employed the no- 
tation rjp = {1 ~ (— 1)P). As we are interested in the 
time-evolution of the initial Fock state, we abbreviated 
unimportant terms as \(j)a) which are the states beside 
the Fock state that are directly coupled via the kinetic 
term to the states \(pp)- 

Using these eigenenergies and states to first order, we 
can now write the time evolving state |'0(^)) &s 



Jy/n{n + 1) 



U 



^(|</.(m,l)) + |0(m,-l))) 



^^^|±lH5:,,sin(.p/L)e-'**|0°). (A5) 



U 



sin(7rp/L)|n) + J/u'Y^l'fa) This formula corresponds to the expression (62| which 



(A3) 



with corresponding energies 

Ep^U - 2{2n + 1) J cos(7rp/i) . (A4) 



one obtains in the unconstrained fermionic approach 
up to first order. As discussed in the unconstrained 
fermionic approach this expression gives a lot of insight 
into the formation and propagation of singly and dou- 
bly occupied sites and can be used to compute all the 
observables that we are interested in. 
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